package parse_bl2seq;

#  get_alignment_data('/share/tools/blast/bl2seq', $seq1, $seq2
#  take the @results[4,5]
#  expand cDNA using expand_cdna_with_alignment_array
#  e.g. expand_cdna_with_alignment_array($cdna_seq, '-', $results[4]);
#  call once prepare_calculate_yn00
#	call in loop				calculate_yn00
#  call when finished					cleanup_calculate_yn00


require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(	get_alignment_data
					expand_cdna_with_alignment_str
					expand_cdna_with_alignment_array
					expand_prot_with_alignment_str
					expand_prot_with_alignment_array
					pos_in_aligned_seq
					aligned_pos_in_orig_seq
					get_alignment_marks
					gap_positions_array_to_str
					gap_positions_str_to_array
					get_best_alignment
					do_get_best_alignment
					get_pi_sequence_acc_code
					BL_BIT_SCORE
					BL_ALIGN_BEG
					BL_ALIGN_LEN
					BL_E___VALUE
					BL_ALN_GAPS1
					BL_ALN_GAPS2
					BL_PI1______
					BL_PI2______
					BL_ENSP1____
					BL_ENSP2____
					BL_TOTAL_LEN
					BL_COVERAGE_
					BL_PROT_SEQ1
					BL_PROT_SEQ2
					BL_CDNA_SEQ1
					BL_CDNA_SEQ2
					YN00_KAKS
					YN00_KA__
					YN00_KS__
					YN00_A___
					YN00_S___
					prepare_calculate_yn00
					calculate_yn00
					cleanup_calculate_yn00
					);

use strict;
use warnings;
use IO::File;
use POSIX qw(tmpnam);
use threads;
use threads::shared;

# function declarations
sub parse_blast_results(\@$$);

# get blast alignment given two gene_ids of specified species and blast options
# $dbh, $gene_id1, $gene_id2, $tax_name1, $tax_name2, $options
# $options = 	{
#						'seg'		=> 1,
#						'coils'		=> 1,
#						'e_value'	=> 1e-4,
#						'blosom'	=> 80,
#						'no_sequences'=> undef,

sub get_best_alignment($$$$);



my $tmp_file_name1;
my $tmp_file_name2;
my $tmp_file_name3;
my $tmp_file_name4;

sub delete_tmp_files
{
	my $tmpfile_mutex :shared;
	if ($tmp_file_name1)
	{
		lock ($tmpfile_mutex);
		unlink $tmp_file_name1 or die "Couldn't unlink $tmp_file_name1 : $!";
		$tmp_file_name1 = undef;
	}
	if ($tmp_file_name2)
	{
		lock ($tmpfile_mutex);
		unlink $tmp_file_name2 or die "Couldn't unlink $tmp_file_name2 : $!";
		$tmp_file_name2 = undef;
	}
	if ($tmp_file_name3)
	{
		lock ($tmpfile_mutex);
		unlink $tmp_file_name3 or die "Couldn't unlink $tmp_file_name3 : $!";
		$tmp_file_name3 = undef;
	}
	if ($tmp_file_name4)
	{
		lock ($tmpfile_mutex);
		unlink $tmp_file_name4 or die "Couldn't unlink $tmp_file_name4 : $!";
		$tmp_file_name4 = undef;
	}

}

#_________________________________________________________________________________________


#	line_wrap

#			puts a newline every n letters
#			n.b. Do not use for constant parameters eg "$3"

#_________________________________________________________________________________________
sub line_wrap($$)
{
	my @array;
	my $seq = $_[0];
	while (length $seq)
	{
		push @array, substr $seq, 0, $_[1], "";
	}
	return join "\n", @array;
}


# delete files automatically when we exit or die
END
{
	delete_tmp_files();
}	


#_________________________________________________________________________________________

#	get_alignment_data

#		input:

#			$blast_prog			name and path of the blast programme
#			$seq1				sequence in fasta format
#			$seq2				sequence in fasta format

#		returns:

		
#			$score				bitscore
#			$alignment_start	start of the alignment
#			$alignment_length	length of the alignment
#			$e-value			e-value
#			$i_pos				reference to an array of position / length pairs
#			$j_pos				reference to an array of position / length pairs
#			$final_len			length of alignment positions + non_aligned positions
#			$coverage			alignment coverage = $alignment_length / $final_len

#_________________________________________________________________________________________
sub get_alignment_data
{
	my ($blast_prog, $seq1, $seq2, $options) = @_;

	die unless length($seq1);
	die unless length($seq2);
	my $seg = ($options->{'seg'} ? '-F t' : '-F f');
	my $blosom = ($options->{'blosom'} ? "-M BLOSUM$options->{'blosom'}" : '');
	my $e_value = ($options->{'e_value'} ? "-e $options->{'e_value'}" : '');
	my $coils = $options->{'coils'};

	my $len1 = length($seq1);
	my $len2 = length($seq2);
	$seq1 = ">seq1\n".line_wrap($seq1,60)."\n";
	$seq2 = ">seq2\n".line_wrap($seq2,60)."\n";
	
	
	#
	#	prepare sequence files for bl2seq
	#
	{
		my $tmpfile_mutex :shared;
		my ($tmp_fh3, $tmp_fh4);
		lock ($tmpfile_mutex);
		do {$tmp_file_name3 = tmpnam()}
			until $tmp_fh3 = IO::File->new($tmp_file_name3, O_RDWR | O_CREAT | O_EXCL );
		do {$tmp_file_name4 = tmpnam()}
			until $tmp_fh4 = IO::File->new($tmp_file_name4, O_RDWR | O_CREAT | O_EXCL );
		if (!$coils)
		{
			print $tmp_fh3 $seq1 or die "Failed to print\n$@";
			print $tmp_fh4 $seq2 or die "Failed to print\n$@";
		}
	}
	
	# open temp files
	if ($coils)
	{
		my ($tmp_fh1, $tmp_fh2);
		{
			my $tmpfile_mutex :shared;
			lock ($tmpfile_mutex);
			do {$tmp_file_name1 = tmpnam()}
				until $tmp_fh1 = IO::File->new($tmp_file_name1, O_RDWR | O_CREAT | O_EXCL );
			do {$tmp_file_name2 = tmpnam()}
				until $tmp_fh2 = IO::File->new($tmp_file_name2, O_RDWR | O_CREAT | O_EXCL );
		}
		print $tmp_fh1 $seq1 or die "Failed to print\n$@";
		print $tmp_fh2 $seq2 or die "Failed to print\n$@";
		close $tmp_fh1;
		close $tmp_fh2;
		system "/share/tools/coils/ncoils -f  < $tmp_file_name2  1> $tmp_file_name4 2> /dev/null";
		system "/share/tools/coils/ncoils -f  < $tmp_file_name1  1> $tmp_file_name3 2> /dev/null";
	}
	
	# blast away
	my @results = split/\n/, `$blast_prog -i $tmp_file_name3 -j $tmp_file_name4 $blosom -X 50 -p blastp $seg $e_value`;
	
	# can remove temp files now
	delete_tmp_files();

	# analyse blast results
	return (parse_blast_results(@results, $len1, $len2));

}

#_________________________________________________________________________________________

#	parse_blast_results

#		input:

#			@results
#				array of lines with the blast results

#		returns:

#			$score				bitscore
#			$alignment_length	length of the alignment
#			$i_pos				reference to an array of position / length pairs
#			$j_pos				reference to an array of position / length pairs

#_________________________________________________________________________________________
sub throw_bad_data($\@)
{
   	die "Error!\nUnexpected data on line ".$_[0]." of the blast output\n"
			.$_[1][$_[0]]."\n";
}

sub gap_positions_array_to_str(\@)
{
	return join ("|", map{$_->[0].':'.$_->[1]} @{$_[0]});
}
sub gap_positions_str_to_array($)
{
	return map{[split (/:/,$_)]} split (/,/, $_[0]);
}
sub parse_blast_results(\@$$)
{
	my ($results, $len_a, $len_b) = @_;
	my $score;
	my ($e_value);
	my ($alignment_length, $alignment_start);
	my (@i_pos, @j_pos);
	for (my $i = 0; $i < @$results; ++$i)
	{
		# get bit score		Score =  703 bits (1814), Expect = 0.0
		next unless ($results->[$i] =~ /^\s*Score\s*=\s*(\d+(?:\.\d+)?).*Expect\s*=\s*([0-9\-\.Ee]+)/);
		$score = $1;
		$e_value = $2;
		$e_value =~ s/^e/1e/;
		
		++$i;
		++$i while (defined $results->[$i] && !length($results->[$i]));
		die "Error!\nPremature end after line $i\n" unless defined $results->[$i];

	
		# get the length of the alignment i.e. "Identities = 527/625 (84%)"
			throw_bad_data($i, @$results) unless
		$results->[$i] =~ /^\s*Identities\s*=\s*\d+\/(\d+)/;
		$alignment_length = $1;
		
		++$i;
		++$i while (defined $results->[$i] && !length($results->[$i]));
		die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
		# go through query/align/subject line triplets of the top alignment
		my $do_initial_align = 1;
		while ($results->[$i] =~ /^\s*Query:/)
		{
			# query
				throw_bad_data($i, @$results) unless
			$results->[$i] =~ /^\s*					# whitespace
								Query:				# name of line
								\s*                 # whitespace
								(\d+)               # start position of query
								\s*                 # whitespace
								([A-Z\-+]+)/x;		# text string
	
			my $i_aa_line	= $2;
			my $i_seq_start = $1;
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	

			#skip past alignment line
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
	
			# subject
			throw_bad_data($i, @$results) unless
				$results->[$i] =~ /^\s*					# whitespace
									Sbjct:				# name of line
									\s*                 # whitespace
									(\d+)               # start position of query
									\s*                 # whitespace
									([A-Z\-+]+)/x;		# text string
			my $j_aa_line = $2;
			my $j_seq_start = $1;
			++$i;
			++$i while (defined $results->[$i] && !length($results->[$i]));
			die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
	
			# note the offset to the first aligned position
			if ($do_initial_align)
			{
	
				$alignment_start = $i_seq_start;
				if ($i_seq_start > $j_seq_start)
				{
					push(@j_pos, [0, $i_seq_start-$j_seq_start]);
				}
				elsif ($i_seq_start < $j_seq_start)
				{
					push(@i_pos, [0, $j_seq_start-$i_seq_start]);
					$alignment_start = $j_seq_start;
				}
				$do_initial_align = 0;
			}
	

			# save the number and position of the aligning characters for each sequence
			while ($i_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($i_aa_line) - $len + $i_seq_start - 1;
				if (@i_pos && $i_pos[-1][0] == $pos)
				{
					$i_pos[-1][1] += $len;
				}
				else
				{
					push(@i_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$i_seq_start -= $len;
			}
	
			while ($j_aa_line =~ /(\-+)/g)
			{
				# save position and length of aligning markers
				my $len = length($1);
				my $pos = pos($j_aa_line) - $len + $j_seq_start - 1;
				if (@j_pos && $j_pos[-1][0] == $pos)
				{
					$j_pos[-1][1] += $len;
				}
				else
				{
					push(@j_pos, [$pos, $len]);
				}
				# adjust starting position to account for the markers already counted
				$j_seq_start -= $len;
			}
		}
		
		if (defined $score)
		{
		#find length of two sequences
			while ($results->[$i] !~ /^length of query:\s*(\d+)/)
			{
				++$i;
				++$i while (defined $results->[$i] && !length($results->[$i]));
				die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
			}
			$len_a = $1;
			while ($results->[$i] !~ /^length of database:\s*(\d+)/)
			{
				++$i;
				++$i while (defined $results->[$i] && !length($results->[$i]));
				die "Error!\nPremature end after line $i\n" unless defined $results->[$i];
			}
			$len_b = $1;
		}
		my $final_len_a = $len_a;
		$final_len_a += $_->[1] for (@i_pos);
		my $final_len_b = $len_b;
		$final_len_b += $_->[1] for (@j_pos);

		# add end alignment marks
		my $final_len = $final_len_a;
		if ($final_len_a > $final_len_b)
		{
			push(@j_pos, [$len_b, $final_len_a - $final_len_b]);
			$final_len = $final_len_a;
		}
		elsif ($final_len_b > $final_len_a)
		{
			push(@i_pos, [$len_a, $final_len_b - $final_len_a]);
			$final_len = $final_len_b;
		}
		return ($score, $alignment_start, $alignment_length, $e_value,
					gap_positions_array_to_str(@i_pos),
					gap_positions_array_to_str(@j_pos),
					$final_len,
					$alignment_length / $final_len);
	}
	return;

}




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	print the modified sequences

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#_________________________________________________________________________________________

#	expand_cDNA_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs in string format 0:63|119:6|263:1

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_cdna_with_alignment_array($$\@)
{
	my ($seq, $char, $align) = @_;
	$char = '-' unless defined $char && length $char;
	$char = substr($char, 0,  1);
	$char = $char x 3;
	
	foreach (reverse @$align)
	{
		substr($seq, $_->[0] * 3, 0) = ($char x $_->[1]);
	}
	return $seq;

}

sub expand_cdna_with_alignment_str($$$)
{
	my ($seq, $char, $align_str) = @_;
	my @align = map {[split '\:', $_]; } split '\|', $align_str;
	return expand_cdna_with_alignment_array($seq, $char, @align)
}


#_________________________________________________________________________________________

#	expand_sequence_to_alignment

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_prot_with_alignment_array($$\@)
{
	my ($seq, $char, $align) = @_;
	$char = '-' unless defined $char && length $char;
	$char = substr($char, 0,  1);
	foreach (reverse @$align)
	{
		substr($seq, $_->[0], 0) = ($char x $_->[1]);
	}
	return $seq;
}
sub expand_prot_with_alignment_str($$$)
{
	my ($seq, $char, $align_str) = @_;
#	return unless $align_str;
	my @align = map {[split '\:', $_]; } split '\|', $align_str;
	return expand_prot_with_alignment_array($seq, $char, @align)
}

#_________________________________________________________________________________________

#	aligned_pos_in_orig_seq

#		input:

#			$pos =  aligned position
#			@list of alignment position/length pairs

#		returns:

#			corresponding position in original alignment
#_________________________________________________________________________________________
sub aligned_pos_in_orig_seq($\@)
{
	my ($pos, $align) = @_;
	return 0 if ($pos < 0);
	my $new_pos = $pos;
	foreach (@$align)
	{
		my $start = $_->[0];
		my $len = $_->[1];
		if ($new_pos >= $start + $len)
		{
			$new_pos -= $len;
			next;
		}
		return ($start - 1) if ($new_pos >= $start);
		return $new_pos;
	}
	return $new_pos;

}
#_________________________________________________________________________________________

#	pos_in_aligned_seq

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub pos_in_aligned_seq($\@)
{
	my ($pos, $align) = @_;
	my $new_pos = $pos;
	foreach (@$align)
	{
		$new_pos += $_->[1] if ($pos >= $_->[0]);
	}
	return $new_pos;

}


#_________________________________________________________________________________________

#	get_alignment_marks

#		input:

#			$ aligned protein sequence with no spaces etc.

#		returns:

#			positions of aligning marks
#_________________________________________________________________________________________
sub get_alignment_marks($)
{
	my ($seq) = @_;
	my @align_pos;
	my $cummulative = 0;
	# save the number and position of the aligning characters for each sequence
	while ($seq =~ /(\-+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		push(@align_pos, [pos($seq) - $len - $cummulative, $len]);
		$cummulative += $len;
	}
	return @align_pos
}


















##########################################################################################

#     get_best_alignment

#        get blast alignment given two gene_ids of specified species and blast options

##########################################################################################
sub get_best_alignment($$$$);
#
# constants
#

# data retrieved for each transcript of each gene
use constant BL_PI__ =>  0;
use constant BL_SEQ_ =>  1;
use constant BL_ENSP =>  2;

# for alignment results
use constant BL_BIT_SCORE =>   0;
use constant BL_ALIGN_BEG =>   1;
use constant BL_ALIGN_LEN =>   2;
use constant BL_E___VALUE =>   3;
use constant BL_ALN_GAPS1 =>   4;
use constant BL_ALN_GAPS2 =>   5;
use constant BL_TOTAL_LEN =>   6;
use constant BL_COVERAGE_ =>   7;
use constant BL_PI1______ =>   8;
use constant BL_PI2______ =>   9;
use constant BL_ENSP1____ =>  10;
use constant BL_ENSP2____ =>  11;
use constant BL_PROT_SEQ1 =>  12;
use constant BL_PROT_SEQ2 =>  13;
use constant BL_CDNA_SEQ1 =>  14;
use constant BL_CDNA_SEQ2 =>  15;



#_________________________________________________________________________________________

#     do_get_best_alignment

#        get the best alignment in all vs all of ortholog transcripts
#_________________________________________________________________________________________
sub do_get_best_alignment($$$)
{
	# arrays of [BL_PI__, BL_SEQ_, BL_ENSP] for each taxon
	my ($all_taxon1_data, $all_taxon2_data, $options) = @_;
	
	my $orig_e_value = $options->{'e_value'};

	#
	#	all against all comparison
	#
	my @best_alignment;
	my $best_combined_len;
	
	for my $taxon1 (@$all_taxon1_data)
	{
		for my $taxon2 (@$all_taxon2_data)
		{
			# ignore if either sequence is too short
			next if (length($taxon1->[BL_SEQ_]) < 50 ||
					 length($taxon2->[BL_SEQ_]) < 50);
			
			# call blast
		#	print STDERR $taxon1->[BL_SEQ_], "\n",
		#				$taxon2->[BL_SEQ_], "\n";
			my @alignment = get_alignment_data("/share/tools/blast/bl2seq",
												$taxon1->[BL_SEQ_],
												$taxon2->[BL_SEQ_],
												$options);

			# make sure blast succeeded before proceeding further
			next unless defined $alignment[BL_BIT_SCORE] &&
						defined $alignment[BL_E___VALUE];
		
			# ignore if e value too high
			next if ($alignment[BL_E___VALUE] > $options->{'e_value'});
			
			# only overwrite current alignment if better bit score,
			# or if equal bit score lower e_value and shorter length
			if (@best_alignment)
			{
				# do not save if score is less
				next if ($alignment[BL_BIT_SCORE] < $best_alignment[BL_BIT_SCORE]);
				
				# if same bitscore
				if ($best_alignment[BL_BIT_SCORE] == $alignment[BL_BIT_SCORE])
				{
					# only save if lower e_value
					next if ($alignment[BL_E___VALUE] > $best_alignment[BL_E___VALUE]);
					
					#  or shorter length
					next if ($alignment[BL_TOTAL_LEN]  < $best_alignment[BL_TOTAL_LEN]);
				}
			}
		
			# use the best e-value as the new cutoff
			$options->{'e_value'} = ($alignment[BL_E___VALUE] == 0 ? 1e-40 :
									 $alignment[BL_E___VALUE]);
			
			#save
			@best_alignment = @alignment;
			$best_alignment[BL_PI1______] = $taxon1->[BL_PI__];
			$best_alignment[BL_PI2______] = $taxon2->[BL_PI__];
			$best_alignment[BL_ENSP1____] = $taxon1->[BL_ENSP];
			$best_alignment[BL_ENSP2____] = $taxon2->[BL_ENSP];
			
			$best_alignment[BL_PROT_SEQ1] = $taxon1->[BL_SEQ_];
			$best_alignment[BL_PROT_SEQ2] = $taxon2->[BL_SEQ_];
		}
	}
	
	# restore original e-value cutoff
	$options->{'e_value'} = $orig_e_value;
	return @best_alignment;
}





#_________________________________________________________________________________________

#     get_cdna_sequence

#_________________________________________________________________________________________
sub get_cdna_sequence($$$)
{
	my ($dbh, $gene_id, $ensp) = @_;
	my $sql_cmd = <<"PL/SQLCMD"; $sql_cmd =~ s/\t/ /;
	SELECT cdna
		FROM cdna_sequences
		NATURAL JOIN ens_id
		WHERE
				gene_id = '$gene_id' AND
				acc_code = '$ensp';
PL/SQLCMD
		
	# get all pi and transcript sequences for the first gene in this ortholog pair
	my @row = $dbh->selectrow_array( $sql_cmd) 				or die DBI::errstr();
	return $row[0];
}




#_________________________________________________________________________________________

#     get_pi_sequence_acc_code

#_________________________________________________________________________________________
sub get_pi_sequence_acc_code($$)
{
	my ($dbh, $gene_id) = @_;
	my $sql_cmd = <<"PL/SQLCMD"; $sql_cmd =~ s/\t/ /;
	SELECT DISTINCT ON (pi) pi, sequence, acc_code
		FROM sequences
		NATURAL join ens_id
		WHERE gene_id = 'GENE_ID';
PL/SQLCMD
		
	# get all pi and transcript sequences for this ortholog pair
	$sql_cmd =~ s/GENE_ID/$gene_id/;
	my $sth = $dbh->prepare( $sql_cmd) 							or die DBI::errstr();
	$sth->execute()												or die DBI::errstr();
	return $sth->fetchall_arrayref([]);
}


#_________________________________________________________________________________________

#     get_best_alignment

#        get blast alignment given two gene_ids of specified species and blast options

#_________________________________________________________________________________________
sub get_best_alignment($$$$)
{
	my ($dbh, $gene_id1, $gene_id2, $options) = @_;

	my ($pi_seq_ensp1, $pi_seq_ensp2);
	if (!defined $options)
	{
		$options = 	{
						'seg'		=> 1,
						'coils'		=> 1,
						'e_value'	=> 1e-4,
						'blosom'	=> 80,
						'no_sequences'=> undef,
					};
	}
	
	# get all pi and transcript sequences for each species of this ortholog pair
	$pi_seq_ensp1 = get_pi_sequence_acc_code($dbh, $gene_id1);
	$pi_seq_ensp2 = get_pi_sequence_acc_code($dbh, $gene_id2);
	my @alignment_data = do_get_best_alignment($pi_seq_ensp1, $pi_seq_ensp2, $options);
	unless (@alignment_data)
	{
		return;
	}
	die unless defined $alignment_data[BL_PI2______];

	if ($options->{no_sequences})
	{
		return (@alignment_data);
	}

	#
	#	return more information about the alignment
	#

	# aligned proteins as well as the aligned region size
	my $prot1 = expand_prot_with_alignment_str( $alignment_data[BL_PROT_SEQ1], '-',
												$alignment_data[BL_ALN_GAPS1]);
	$alignment_data[BL_PROT_SEQ1] = substr($prot1, $alignment_data[BL_ALIGN_BEG],
												$alignment_data[BL_ALIGN_LEN]);
	my $prot2 = expand_prot_with_alignment_str( $alignment_data[BL_PROT_SEQ2], '-',
												$alignment_data[BL_ALN_GAPS2]);
	$alignment_data[BL_PROT_SEQ2] = substr($prot2, $alignment_data[BL_ALIGN_BEG],
												$alignment_data[BL_ALIGN_LEN]);
	
	
	# cDNA
	$alignment_data[BL_CDNA_SEQ1] = substr(expand_cdna_with_alignment_str(
											get_cdna_sequence($dbh, $gene_id1,
															  $alignment_data[BL_ENSP1____]),
											'-',
											$alignment_data[BL_ALN_GAPS1]),
										$alignment_data[BL_ALIGN_BEG] * 3,
										$alignment_data[BL_ALIGN_LEN] * 3);

	$alignment_data[BL_CDNA_SEQ2] = substr(expand_cdna_with_alignment_str(
											get_cdna_sequence($dbh, $gene_id2,
															  $alignment_data[BL_ENSP2____]),
											'-',
											$alignment_data[BL_ALN_GAPS2]),
										$alignment_data[BL_ALIGN_BEG] * 3,
										$alignment_data[BL_ALIGN_LEN] * 3);
	return (@alignment_data);

}

sub prepare_calculate_yn00($)
{
	my ($dir) = @_;
	mkdir ("$dir/yn00_tmp");
	open FH_YN00_CTL, ">$dir/yn00_tmp/yn00.ctl";
	print FH_YN00_CTL <<FILE;
seqfile = seq.tmp
outfile = yn
verbose = 0
noisy = 0
icode = 0
weighting = 1
commonf3x4 = 0
FILE
	close FH_YN00_CTL;
}

sub cleanup_calculate_yn00($)
{
	my ($dir) = @_;
	unlink(glob("$dir/yn00_tmp/*"));
	rmdir("$dir/yn00_tmp");
}

# for alignment results
use constant YN00_KAKS =>  0;
use constant YN00_KA__ =>  1;
use constant YN00_KS__ =>  2;
use constant YN00_A___ =>  3;
use constant YN00_S___ =>  4;

sub calculate_yn00($$$)
{
	my ($dir, $seq1, $seq2) = @_;
	die unless length $seq1 == length $seq2;
	open FH_SEQ, ">$dir/yn00_tmp/seq.tmp";
	print FH_SEQ "2 ".length($seq1), "\n";
	print FH_SEQ "a\n", $seq1, "\n";
	print FH_SEQ "b\n", $seq2, "\n";
	close FH_SEQ;
	unlink "$dir/yn00_tmp/yn";
	
	chdir "$dir/yn00_tmp";
	system ("/share/tools/paml/yn00 < yn00.ctl > /dev/null");
	chdir "..";
	unless (open FH_YN, "$dir/yn00_tmp/yn")
	{
		return;
	}
	while (<FH_YN>)
	{
		next unless /seq\.\s+seq./;
		<FH_YN>;
		$_ = <FH_YN>;
		my @values = split /[\s\+\-]+/;
		shift @values;
		return unless @values == 11;
		my @results = ($values[6],                 # kaks
					   $values[7],                 # ka
					   $values[9],                 # Ks
					   $values[3],                 # A
					   $values[2]                  # S
					   );
	#	print STDERR join "\n", @results, "\n";
		# ignore invalid values
		for (@results)
		{
			unless (defined $_ && /\d/ && $_ >= 0.0)
			{
				$_ = '\N' ;
			}
		}
		
		# ignore out of range
		for (YN00_KAKS, YN00_KA__, YN00_KS__)
		{
			$results[$_] = '\N'  unless ($results[$_] < 50.0);
		}
		return (@results);				
	}
	return;
}





1;

